xxxxxxxxxxIn this notebook, we're solving the simple pendulum problem analytically and numerically. The simple pendulum is a classical example of a nonlinear, second-order, ordinary differential equation (ODE).In this notebook, we're solving the simple pendulum problem analytically and numerically. The simple pendulum is a classical example of a nonlinear, second-order, ordinary differential equation (ODE).
xxxxxxxxxx## Problem StatementThe simple pendulum is a mass $m$ attached to a massless, frictionless rod of length $l$ that is free to rotate about a fixed point. The pendulum is released from rest at an angle $\theta_0$ from the verticalThe simple pendulum is a mass attached to a massless, frictionless rod of length that is free to rotate about a fixed point. The pendulum is released from rest at an angle from the vertical
xxxxxxxxxx## Necessary libraries* `numpy` - for numerical calculations* `matplotlib` - for plotting* `scipy.integrate` - for solving ODEs* `sympy` - for symbolic calculationsnumpy - for numerical calculationsmatplotlib - for plottingscipy.integrate - for solving ODEssympy - for symbolic calculationsxxxxxxxxxxRun the next cell in case you don't have these libraries installed.Run the next cell in case you don't have these libraries installed.
%pip install sympy numpy matplotlib scipy --quietimport numpy as npimport matplotlib.pyplot as pltfrom scipy.integrate import odeintfrom matplotlib import animationfrom IPython.display import HTMLimport sympy as symxxxxxxxxxxWe need to define the following constants and variables:* $g$ - acceleration due to gravity* $l$ - length of the pendulum* $m$ - mass of the pendulum* $\theta$ - angle of the pendulum from the vertical (function of time)* $t$ - timeWe need to define the following constants and variables:
g = sym.symbols('g')l = sym.symbols('l')m = sym.symbols('m')t = sym.symbols('t')theta = sym.symbols(r'\theta', cls=sym.Function)theta = theta(t)xxxxxxxxxxNow we need to define $x$ and $y$ in terms of $\theta$.* $x$ - horizontal position of the pendulum - $x = l \sin(\theta)$* $y$ - vertical position of the pendulum - $y = -l \cos(\theta)$Now we need to define and in terms of .
x = l * sym.sin(theta)y = -l * sym.cos(theta)xxxxxxxxxxOnce we have $x$ and $y$, we can get the kinetic and potential energies of the pendulum.* $T$ - kinetic energy of the pendulum $$T = \frac{1}{2} m \dot{x}^2 + \frac{1}{2} m \dot{y}^2$$* $V$ - potential energy of the pendulum $$V = m g y$$Once we have and , we can get the kinetic and potential energies of the pendulum.
# Kinetic energyT = sym.Rational(1, 2) * m * (x.diff(t) ** 2 + y.diff(t) ** 2)# Potential energyV = m * g * yxxxxxxxxxxNow, we can get the Lagrangian of the pendulum.* $L$ - Lagrangian of the pendulum - $L = T - V$Now, we can get the Lagrangian of the pendulum.
L = T - VL.simplify()xxxxxxxxxxWith the Lagrangian, we can get the equations of motion.$$\frac{d}{dt} \frac{\partial L}{\partial \dot{\theta}} - \frac{\partial L}{\partial \theta} = 0$$With the Lagrangian, we can get the equations of motion.
LE = (sym.diff(sym.diff(L, sym.diff(theta, t)), t) - sym.diff(L, theta)).simplify()xxxxxxxxxxNow, let's solve the equation for $\ddot{\theta}$.Now, let's solve the equation for .
sol = sym.solve(LE, sym.diff(theta, t, t), simplify=False, rational=False)sol[0]xxxxxxxxxxWith the analitycal solution for $\ddot{\theta}$, we can get the numerical solution for $\dot{\theta}$ and $\theta$.With the analitycal solution for , we can get the numerical solution for and .
sol_f = sym.lambdify((t, theta, sym.diff(theta, t), g, l, m), sol[0])xxxxxxxxxxLet's define a vector $\vec{S}$ that returns $\frac{d\vec{S}}{dt}$.Let's define a vector that returns .
def dSdt(S, t, g, l, m): theta, omega = S dtheta_dt = omega domega_dt = sol_f(t, theta, omega, g, l, m) return dtheta_dt, domega_dtxxxxxxxxxxNow, let's plot the results.Now, let's plot the results.
t = np.linspace(0, 40, 1001)theta0 = np.pi * 0.12theta_dot0 = 0.5g = 9.81l = 1m = 1ans = odeint(dSdt, (theta0, theta_dot0), t, args=(g, l, m))# Plot of the position and velocity of the pendulumfig, ax = plt.subplots(2, 1, figsize=(10, 10))ax[0].plot(t, ans[:, 0])ax[0].set_xlabel('Time (s)')ax[0].set_ylabel('Angle (rad)')ax[1].plot(t, ans[:, 1])ax[1].set_xlabel('Time (s)')ax[1].set_ylabel('Angular velocity (rad/s)')plt.show()xxxxxxxxxxPlot phase space using streamplots.Plot phase space using streamplots.
fig, ax = plt.subplots(1, 1, figsize=(10, 5))x_values, y_values = np.meshgrid(np.linspace(-10, 10, 20), np.linspace(-10, 10, 20))u_values = y_valuesv_values = -g/l * np.sin(x_values)ax.streamplot(x_values, y_values, u_values, v_values)ax.set_xlabel('Angle (rad)')ax.set_ylabel('Angular velocity (rad/s)')ax.grid()plt.show()xxxxxxxxxxTo create an animation, we need to define a function that returns $x$ and $y$ coordinates for each frame.To create an animation, we need to define a function that returns and coordinates for each frame.
def get_x_y(theta, l): x = l * np.sin(theta) y = -l * np.cos(theta) return x, y# Animation of the pendulumfig, ax = plt.subplots(figsize=(5, 5))ax.set_xlim(-1.1, 1.1)ax.set_ylim(-1.1, 1.1)ax.set_aspect('equal')ax.grid()line, = ax.plot([], [], 'o-', lw=2)time_template = 'time = %.1fs'time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)# X and Y coordinates of the pendulumx, y = get_x_y(ans[:, 0], l)def init(): line.set_data([], []) time_text.set_text('') return line, time_textdef animate(i): thisx = [0, x[i]] thisy = [0, y[i]] line.set_data(thisx, thisy) time_text.set_text(time_template % (i * t[1])) return line, time_textani = animation.FuncAnimation(fig, animate, np.arange(1, len(ans)), interval=25, blit=True, init_func=init)plt.rcParams['animation.embed_limit'] = 100HTML(ani.to_jshtml())